Control of particle clustering in turbulence by polymer additives 
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We study the clustering properties of inertial particles in a turbulent viscoelastic fluid. The in- 
vestigation is carried out by means of direct numerical simulations of turbulence in the Oldroyd-B 
model. The effects of polymers on the small scale properties of homogeneous turbulence are con- 
sidered in relation with their consequences on clustering of particles, both lighter and heavier than 
the carrying fluid. We show that, depending on particle and flow parameters, polymers can either 
increase or decrease clustering. 



Clustering of inertial particles in turbulent flows is rel- 
evant for meteorology and engineering, as well as fun- 
damental research. It is believed to play a crucial role 
in rain-drop formation [l|, as well as in the aggregation 
of proto-planetesimals in Keplerian accretion disks [2j. 
The physical mechanism which originates such cluster- 
ing is indeed rather simple: particles heavier than the 
fluid in which they are transported experience inertial 
forces which expel them from vortices; particles lighter 
than the fluid are attracted into vortical structures, for 
similar reasons 04a|- In realistic flows, however, parti- 
cles are advected by the small scale vortical structures of 
turbulent flows: these have highly non-trivial statistical 
features, resulting in a complex clustering process which 
is still far from being completely understood. From the 
point of view of applications, the properties of concentra- 
tion and distribution of inertial particles play a crucial 
role in engineering and for the design of industrial pro- 
cesses involving combustion and mixing [6H&]. Suspen- 
sions of particles in viscoelastic fluids are used in many 
products of commercial and industrial relevance [9( . 

In this paper we investigate, by means of direct numer- 
ical simulations of a turbulent flow, how the clustering 
properties of a dilute suspension of inertial particles can 
be affected by the addition of small amounts of polymer 
additives. The effects induced by polymers on turbulent 
flows are themselves of enormous relevance. It is enough 
to mention the celebrated drag reduction effect which 
occurs in pipe flows lldl . or the recently discovered elas- 
tic turbulence regime [lj| . Polymers have striking effects 
also on Lagrangian properties of the flow. In particular it 
has been shown that polymer addition in turbulent flows 
reduces the chaoticity of Lagrangian trajectories [12| and 
affects acceleration of fluid tracers [l3| . Conversely in the 
elastic turbulence regime polymers are able to generate 
Lagrangian chaos in flows at vanishing Reynolds num- 
ber, which would be non chaotic in the Newtonian case 

nam. 

Here we show that the addition of polymers in a tur- 
bulent flow has important effects on the statistical prop- 
erties of inertial particles which can result in both an in- 
crease or a decrease of the clustering. An example of the 
effect of polymers on clustering is shown in Fig. [1] which 
represents the distribution of an ensemble of inertial par- 



ticles in a turbulent flow before and after the introduction 
of polymers. It is evident, already at the qualitative level 
of Fig.[IJ that polymers are able to change the statistical 
distribution of particles. We show that these effects can 
be understood and quantified in terms of the Lyapunov 
exponents of inertial particles, which are very sensitive to 
the presence of polymers. Previous systematic investiga- 
tions of inertial particle dynamics in Newtonian turbulent 
flows [15| and stochastic flows [16| have shown that clus- 
tering (quantified by means of the Lyapunov Dimension 
of particle attractor) is maximum when the particle re- 
laxation time is of the order of the shortest characteristic 
time of the flow. 

We consider the case of a dilute suspension of small 
inertial particles, in which the effects of the disturbance 
flow induced by the particles can be neglected. The dy- 
namics of the suspension is hence modeled by an ensem- 
ble of non-interacting point particles, which experience 
viscous drag and added mass forces. The equation of 
motion of each particle reads [18| : 



dx 

It 
dv 

~dl 



-L[v-u(x(t),t)]+p^ 
t s at 



(1) 
(2) 



where t$ — a /(3/3f) is the Stokes relaxation time, a is 
the particle radius, j3 = 3pf/(pf + 2p p ) (p p and pf repre- 
senting particle and fluid densities respectively) and v is 
the kinematic viscosity of the fluid (replaced by the to- 
tal viscosity vt in a viscoelastic fluid, see below). Light 
(heavy) particles correspond to /3 > 1 (/3 < 1). In this 
work we consider the two extreme cases of very light par- 
ticles (e.g. air bubble in water) for which /3 — 3 and very 
heavy particles with j3 — 0. We define the Stokes number 
as St = TgAj, where A° is the maximum Lyapunov expo- 
nent of neutral Lagrangian tracers (i.e. St = particles) 
in the flow. With this definition, maximum clustering is 
obtained for St ~ 0.1 [H, [H|. 

The viscoelastic flow u(x, t) in which the particles 
are suspended can be described by standard viscoelas- 
tic models, such as the Oldroyd-B model or the nonlin- 
ear FENE-P model, which accounts for the finite exten- 
sibility of polymers. In spite of their simplicity, these 
models are able to reproduce many relevant properties of 
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FIG. 1: Section on plane z — of the distribution of heavy 
particles with rg = 0.035 (upper panels) and light particles 
with rs — 0.03 (lower panels) in statistically stationary con- 
ditions in a Newtonian flow (left) and a viscoelastic flow at 
Wi = 1 (right). Both flows are forced with the same forc- 
ing f (x, t) 5-correlated in time and localized on large scales. 
Numerical simulations are done by a pseudo-spectral, fully 
dealiased code at resolution 256 . For the viscoelastic sim- 
ulations, a small diffusive term is added to Q to prevent 
numerical instabilities [17t ] . 



dilute polymer solutions, including turbulent drag reduc- 
tion [Jjl, 120] and elastic turbulence phenomenology [21| . 
Here we choose the Oldroyd-B model (22|, in which the 
coupled dynamics of the velocity field u(x, t) and the 
polymer conformation tensor tx(x, t) (which is propor- 
tional to local square polymer elongation) reads: 

■^Iu-Vm = -Vp + z/V 2 w+ — V -cr + f (3) 
at t v 



d<J 
at 



(Vu) 1 -cr + a- (Vu) (a - I)(4) 



The total viscosity of the solution i/j- = v{l + 7) is writ- 
ten in terms of the kinematic viscosity of the solvent v 
and the zero-shear contribution of the polymer 7 which 
is proportional to the polymer concentration. The poly- 
mer time t p represents the longest relaxation time to the 
equilibrium configuration (a — I in dimensionless units). 
Viscoelasticity of the turbulent flow is parametrized by 
the Weissenberg number Wi, the ratio between t p and a 
characteristic time of the flow. Here we use Wi = TpAf 
where A^ is the Lagrangian Lyapunov exponent of the 
Newtonian flow, before the addition of polymers (i.e. ([3]) 
with 7 = 0). We stress that \\ introduced above refers 
instead to the specific flow that carries the suspension 
and it clearly depends on Wi. Therefore A^ 



In the following we discuss results obtained by in- 
tegrating numerically the viscoelastic model dSJU]) at 
high resolution for different values of Wi (see Table [TJ) . 
The flow is sustained by a stochastic Gaussian forcing 
f (x, t) (5-correlatcd in time and localized on large scales. 
Fluid equations were integrated by means of a standard, 
fully dealiased, pseudo spectral code, on a cubic, triple- 
periodic domain with 256 grid points per side. When the 
flow reaches a turbulent, statistically stationary state, 
different families (i.e. with different values of parameters 
/3 and rg) of inertial particles arc injected, with initial 
homogeneous distribution in space, and their motion in- 
tegrated according to (Q][2]). For each value of Wi, we in- 
tegrated the motion of 1024 particles for each of 21 values 
of rs and two values of /?, namely very heavy particles 
with (3 = and "bubbles" with (3 = 3. 

As an effect of inertia the distribution of particles does 
not remain homogeneous and evolves to a fractal set dy- 
namically evolving with the flow, such as the examples 
shown in Fig. [TJ In the language of dynamical systems, 
the equations (JJ5]) for particle motion represent a dissi- 
pative system whose chaotic trajectories evolve to a frac- 
tal attractor (which evolves in time following the flow). 
A quantitative measure of clustering at small scales is 
therefore obtained by measuring the fractal dimension of 
the attractor (for each family of particles) using the Lya- 
punov dimension [161 123| defined in terms of Lyapunov 
exponents as Dl = K + J2i=i -WI'W+il where K is the 
largest integer for which J2i=i^i — UM- Since the 
space distribution of the particles is the projection of the 
attractor on the sub-space of particle positions, the frac- 
tal dimension of clusters is given by min(Dj,, 3) |25l [26J|. 
provided that the projection is generic (for a discussion 
on this issue see e.g. [27|), This implies that Dl < 3 gives 
fractal distributions of dimension Dl, while Dl > 3 cor- 
responds to space-filling configurations, which however 
can be non-homogeneous. 

In Fig. [5] we plot the fractal dimensions for both heavy 
and light particles as a function of rg for the three simu- 
lations at different Wi. It is evident that the addition of 
polymer changes substantially the clustering properties 
of the particles, both increasing Dl and reducing Dl 
depending on value of t$ . Figure [TJ shows examples of 
clustering reduction, for heavy and light particles respec- 
tively. The upper panels refer to heavy particles (/? = 0) 



Wi 



£/ 



Ai> 



0.28 0.28 0.76 1.36 
0.5 0.28 0.18 0.73 1.08 

1 0.28 0.092 0.68 0.75 



Wwi 



TABLE I: Parameters for the Newtonian and viscoelastic sim- 
ulations. The Weissenberg number Wi, energy input £/, vis- 
cous dissipation rate e„, rms velocity u rma and Lagrangian 
Lyapunov exponent A? of the carrier flow are shown. In both 
viscoelastic runs an additional dissipative term was added on 
polymers (see text), with coefficient v v = 2.3 x 10 3 
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FIG. 3: Energy spectra for the Newtonian case Wi = 
(squares) and for the viscoelastic ones Wi = 0.5 (circles) and 
Wi = 1 (triangles). The depletion due to polymer feed-back 
is evident on large wavenumbers, while the larger scales are 
unaffected. The effect of polymers extends at lower wave- 
numbers as Wi increases. Inset: viscous energy dissipation e„ 
during a typical time interval in the stationary simulations, 
for the Newtonian (solid line), Wi = 0.5 (dashed line) and 
Wi = 1 (dash-dot) flows. The decrease in e v with Wi is 
evident, as well as the reduction in fluctuations. 



FIG. 2: Lyapunov dimension for light (upper panel) and 
heavy (lower panel) particles plotted as a function of ts. Dif- 
ferent lines correspond to the different Weissenberg numbers: 
Wi = (squares), Wi = 0.5 (circles) and Wi = 1.0 (triangles). 



with ts = 0.035, while the bottom ones are extracted 
from a simulation with (3 = 3 and ts = 0.03. Both val- 
ues of Stokes time are, for the Newtonian flow, on the left 
of the minimum in Dl. As a consequence, polymers pro- 
duce a reduction of clustering. Such effect is more visible 
for light particles. A possible reason for this difference 
will be discussed further on. 

The mechanism at the basis of this effect is not triv- 
ial and is a consequence of the change induced by the 
polymers on the small-scale properties of the turbulent 
flow. In Figj3]we plot the energy spectra for the differ- 
ent Wi numbers. The effect of polymers is evident in the 
high-wavenumber range where velocity fluctuations are 
clearly suppressed, resulting in a depletion of the energy 
spectrum, while large-scale fluctuations are unaffected. 

Indeed one can expect that only the fastest eddies 
of the flow, i.e. those whose eddy turn-over time ti is 
shorter that the polymer relaxation time t p , can produce 
a significant elongation of polymers. The elastic feedback 
therefore affects only small scales £ with ti < t p . Con- 
versely, large scales exhibit the same phenomenology of 
a Newtonian flow, characterized by a turbulent cascade 
with a constant energy flux equal to the energy input rate 
£f . The turbulent cascade proceeds almost unaffected by 



the presence of polymers down to the Lumley scale £l, 
whose eddy turn-over time equals the polymer relaxation 
time. A dimensional estimate, based on the Kolmogorov 
scaling for the typical velocity ug ~ eJ i 1 ' 3 and turn- 
over time ti = Ijug ~ el' £ 2 ^ 3 of an eddy of size £, 

3/2 1 /2 

gives £l = t p ' e/ . Polymers would therefore affect 



-/ 



only the small scales £ < £^. Our results are in qualita- 
tive agreement with this picture: the Wi = 0.5 spectrum 
differs from the Newtonian one only for k > 8, while at 
Wi = 1 polymers are active over a larger range of scales. 
The reduction of kinetic energy at small scales, due to 
the transfer of energy to the polymers, is accompanied 
by a reduction of the viscous dissipation e v = v((Vu) 2 ) 
at fixed energy input ef, as can be seen from Table U and 
in the inset of FigO This phenomenon has been previ- 
ously observed both in forced and decaying simulations of 
statistically homogeneous and isotropic turbulence (see, 

e.g., Hi S3). 

The suppression of small-scale motions caused by poly- 
mers has major consequences also on the Lagrangian 
statistics. It is responsible of the reduction of chaoticity 
of Lagrangian trajectories [3fJ. Indeed the chaoticity of 
the flow is directly related to its stretching efficiency via 
the Lyapunov exponents. When polymers are stretched, 
the elastic stress tensor produces a negative feed-back 
on small scale stretching, thus reducing the degree of 
chaoticity of the flow [3(J, [31| . This effect is clearly ob- 
servable in the decrease of the Lagrangian Lyapunov ex- 
ponent of the flow at increasing polymer elasticity (see 
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FIG. 4: Comparison between the Cramer functions of the 
stretching rate 71 computed at Wi = (solid line),Wi = 
0.5 (dashed line),Wi = 1 (dash-dot). Inset: first Lagrangian 
Lyapunov exponent A? (circles) and width /i (squares) of the 
Cramer function (see text) as a function of Wi. The Lyapunov 
exponents are compared with the Newtonian value Ai . 



the inset of FigEJ). 

It is worth to notice that, because of polymers counter- 
action, the Lyapunov exponent of the resulting viscoelas- 
tic flow is smaller than r" 1 . In other words, the Wi num- 
ber computed a posteriori (i.e. after polymer injection) 
is always smaller than unity. This is not in contrast with 
the hypothesis that polymers have a strong active effect 
on the flow mainly when they are stretched, i.e. above 
the so-called coil-stretch transition, which is expected to 
happen around Wi ~ 1 [32]. Indeed, the Lyapunov expo- 
nent simply provides a measure of the average stretching 
in a chaotic flow. One should bear in mind that large 
fluctuations of the stretching rates (and therefore strong 
viscoelastic effects) can occur also when Wi < 1. 

Detailed information on the fluctuations of the stretch- 
ing rates can be obtained from the statistics of the Fi- 
nite Time Lyapunov Exponents (FTLE) 7i . The FTLE 
are defined via the exponential growth rate during a fi- 
nite time T of an infinitesimal M-dimensional volume 
(l/T)ln[V M (T)/V M (0)} 0j. From the 
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definition of the Lyapunov exponents it follows that 
lim/r^oo if = Aj. A large deviation approach suggests 
that the probability density function (PDF) of the largest 
stretching rate 71 measured over a long time T 3> 1/Ai 
takes the asymptotic form Pr (71) ~ N(t) exp[— H(ji)T] 
where the Cramer function H(pf\) is convex and obeys 
the conditions -ff(Ai) = 0, -ff'(Ai) = 0. We computed 
the Cramer function for the Lagrangian FTLE for the 
Newtonian case and the two viscoelastic cases. In the in- 
set of Fig. [4] we plotted the average of the stretching rates 
(i.e., the first Lagrangian Lyapunov exponent of the flow 
A°) and the rescaled variance /i = T(^f), for the three 
values of Wi that we considered. The decrease of the 
Lyapunov exponent (rescaled with the Newtonian value 
Af for comparison) gives a measure of the decrease in 
the chaoticity of the flow, due to the action of Polymers. 
On the other hand, we also observe a decrease in the rela- 
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FIG. 5: Lyapunov dimension for light (upper panel) 

and heavy (lower panel) particles plotted as a function of 
St = rsA?. Different lines correspond to the different Weis- 
senberg numbers with symbols as in Fig. [2] 



tive variance /i/Aj 1 , which implies that polymer feedback 
induces also a reduction of the fluctuations of of stretch- 
ing rates. Inspection of the main panel of Fig. 0] how- 
ever, shows that fluctuations are not reduced uniformly. 
Indeed, the shape of P(7i) changes when polymers are 
added. As is evident in Fig [4] elasticity has the effect of 
raising the right branch of the Cramer function, while the 
left one is comparatively less affected. Given the defini- 
tion of #(71), this amounts to a relative suppression of 
positive fluctuations in the stretching rate: as one could 
expect, polymers have a larger (negative) feedback on 
events of larger stretching. 

The effect of polymers on Lyapunov exponents and 
the Lagrangian nature of the latter suggests to introduce 
the dimensionless Stokes number defined as St = TgA® 
which depends on Wi by the dependence of A? shown in 
Fig. [4j Figure [5] shows the Lyapunov dimension Dl for 
both heavy and light particles as a function of St. It is 
evident that, with respect to Fig. [2J the collapse of the 
curves at different Wi is improved. In particular, the 
minimum of the fractal dimension (which corresponds 
to maximum clustering) occurs almost for the same St 
number. Still, some differences are observable, in par- 
ticular for small St in the case of light particles. This 



can be understood by the following argument. Bubbles, 
at variance with heavy particles, have tendency to con- 
centrate on filaments of high vorticity. Indeed, while the 
minimal dimension for heavy particles is about 2.5 (at 
St ~ 0.1), for light particles at maximal clustering it 
becomes as small as 1.26. Vortex filaments correspond 
to quasi-one-dimensional regions of intense stretching, in 
the direction longitudinal to the vortex, which give ma- 
jor contributions to the right tail of the Cramer function. 
As shown in Fig. [H the effects of polymers on the dis- 
tribution of Lyapunov exponent is more evident in this 
region of strong fluctuations, where the distribution does 
not rescale with A®. It is therefore not surprising that 
also the effects on clustering of light particles cannot be 
completely absorbed in the rescaling of ts with the mean 
stretching rate A J. 

As the fractal dimension is given by a combination of 
the Lyapunov exponents, in order to better understand 
the differences on light and heavy particles, in Fig(B] we 
show the first three Lyapunov exponents as a function 
of St. The first observation is that bubbles, at variance 
with heavy particles, exhibit negative values of A2, con- 
sistently with the lower value of Dt, and the tendency of 
light particles to concentrate towards vortex filaments. 

The first Lyapunov exponent decreases with Wi for any 
value of St, thus indicating that the phenomenon of chaos 
reduction, already discussed for the case of Lagrangian 
tracers, is generic also for inertial particles. On the con- 
trary, the second Lyapunov exponent shows a different 
behavior for light and heavy particles: it increases for 
the former but slightly decreases for the latter. Figure [5] 
shows that the effect of polymers is not a simple rescal- 
ing of the Lyapunov spectrum, which would trivially keep 
the dimension Dl unchanged. From this point of view, 
the almost perfect rescaling of the Lyapunov dimensions 
shown in Fig. [5] is quite surprising and arises as the result 
of compensations of different effects. 

In conclusion, we investigated the clustering proper- 
ties of inertial (heavy and light) particles in a turbulent 
viscoelastic fluid. The main effect of polymers on turbu- 
lent flows is to counteract small-scale fluctuations and to 
reduce its chaoticity. Quantitatively, this results in a de- 
crease in the first Lyapunov exponent of the flow, which, 
in turn, affects clustering of inertial particles. The latter 
can be quantified by means of the fractal (Lyapunov) di- 
mension of particle distributions. Although the effects of 
polymers on the particle Lyapunov exponents are com- 
plex and qualitatively different for light and heavy par- 
ticles, the overall effect on fractal dimension is relatively 



simple and can be rephrased in the rescaling of the char- 
acteristic time of the flow. Indeed, when particle inertia 
is parametrized by the Stokes number St defined with 
the Lyapunov time of the flow, one can approximately 
rescale the curves D^iSt) at all Wi. In contrast, as poly- 
mers do not affect large scale properties of the flow, a 
parametrization of particle inertia based on integral time 
scales would not show a collapse of the curves Dl (St) at 
different Wi. As a consequence, any prediction of par- 
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FIG. 6: The first three Lyapunov exponents for light (/3 = 3, 
upper panel) and heavy (/? = 0, lower panel) particles, at 
different Wi. Continuous, dashed and dotted lines represent 
the first, second and third Lyapunov exponents, while symbols 
correspond to different Wi as in Fig. [2] 



ticle clustering in turbulent polymeric solutions requires 
an accurate estimate of small scale stretching rates. 
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